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■ Abstract — The existence and stability of dissipative breathers in rf SQUID (Superconducting 

^ ' Quantum Interference Device) arrays is investigated numerically. In such arrays, the nonlincarity 

CL|, which is intrinsic to each SQUID, along with the weak magnetic coupling of each SQUID to 

^ ' its nearest neighbors, result in the formation of discrete breathers. We analyze several discrete 

I breather excitations in rf SQUID arrays driven by alternating flux sources in the presence of losses. 

00 . The delicate balance between internal power losses and input power, results in the formation 

of dissipative discrete breather (DDB) structures up to relatively large coupling parameters. 
It is shown that DDBs may locally alter the magnetic response of an rf SQUID array from 
paramagnetic to diamagnetic or vice versa. 
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1. INTRODUCTION 



The discrete breathers (DBs), which are also known as intrinsic localized modes (ILMs), belong to 
a class of nonlinear excitations that appear generically in discrete and spatially extended systems 
[T]. They are loosely defined as spatially localized, time-periodic and stable excitations, that can be 
J> ' produced spontaneously in a nonlinear lattice of weakly coupled elements as a result of fluctuations 

O ! [2], disorded [3|, or by purely deterministic mechanisms @]. The last two decades, a large number 

' of theoretical and experimental studies have explored the existence and the properties of DBs 

in a variety of nonlinear discrete systems. Nowadays, there are rigorous mathematical proofs of 
existence of DBs both for energy conserved and dissipative systems [5l |6], and several numerical 
algorithms for their accurate construction have been proposed [71 18]. Moreover, they have been 
observed experimentally in a variety of systems, including solid state mixed- valence transition metal 
I complexes [9], quasi-one dimensional antiferromagnetic chains [10], arrays of Josephson junctions 

micromechanical oscillators [12J, optical waveguide systems |13j . layered crystal insulator at 
300i^ [14J, and proteins [IS]. 

From the perspective of applications to experimental situations where an excitation is subjected 
to dissipation and external driving, dissipative DBs (DDBs) are more relevant than their energy 
conserved counterparts. The dynamics of DDBs is governed by a delicate balance between the input 
power and internal power losses. Recently, DDBs have been demonstrated numerically in discrete 
and nonlinear magnetic metamaterial (MM) models [16\ [T7] . The MMs are artificial composites that 
exhibit electromagnetic (EM) properties not available in naturally occuring materials. They are 
typically made of subwavelength resonant elements like, for example, the split-ring resonator (SRR). 
When driven by an alternating EM field, the MMs exhibit large magnetic response, either positive 
or negative, at frequencies ranging from the microwave up to the Terahertz and the optical bands 
[191 120j . The magnetic response of materials at those frequencies is particularly important for the 
implementation of devices such as compact cavities, tunable mirrors, isolators, and converters. The 
nonlinearity offers the possibility to achieve dynamic control over the response of a metamaterial in 
real time, and thus tuning its properties by changing the intensity of the external field. Recently, 
the construction of nonlinear SRR-based MMs [18j gives the opportunity to test experimentally the 
existence of DDBs in those materials. 

It has been suggested that periodic rf SQUID arrays can operate as nonlinear MMs in mi- 
crowaves, due to the resonant nature of the SQUID itself and the nonlinearity that is inherent to it 
|21| . The combined effects of nonlinearity and discreteness (also inherent in rf SQUID arrays), may 
lead in the generation of nonlinear excitations in the form of DDBs [22]. In the present work we 
investigate numerically the existence and stability of DDBs in rf SQUID arrays. In the next section 
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Figure 1: Left panel: Schematic drawing of a ring-shaped rf SQUID. Middle panel: Equivalent circuit for an 
rf SQUID in an alternating magnetic field. Right panel: Schematic drawing of a two-dimensional orthogonal 
array of identical rf SQUIDs. 



we shortly describe rf SQUID array model, which consists a simple realization of a planar MM. In 
section 3 we present several types of DDBs that have been constructed using standard numerical 
algorithms, and we discuss their magnetic response. We finish in section 4 with the conclusions. 

2. RF SQUID METAMATERIAL MODEL 

An rf SQUID, shown schematically in the left panel of Fig. 1, consists of a superconducting ring 
interrupted by a Josephson junction (JJ) [23]. When driven by an alternating magnetic field, 
the induced supercurrents in the ring are determined by the J J through the Josephson relations. 
Adopting the resistively and capacitively shunted junction (RCSJ) model for the J J [23], an rf 
SQUID in an alternating field H^xt = H perpendicular to its plane is equivalent to the lumped 
circuit model shown in the middle panel of Fig. 1. That circuit consists of an inductance L in 
series with an ideal Josephson element Ic (i.e., for which I = Jcsin^, where Ic is the critical current 
of the JJ and (p is the Josephson phase) shunted by a capacitor C and a resistor R, driven by an 
alternating flux ^(,xt{H). 

Consider a planar rf SQUID array consisting of identical units (right panel of Fig. 1), arranged 
in an orthogonal lattice with constants and dy in the x and y directions, respectively. That 
system is placed in a uniform magnetic field H = Hue + Hac sin(a;t), where uj is the frequency 
and t is the temporal variable, perpendicular to the SQUID rings. The field induces a supercurrent 
Inm in the nm— th SQUID through the flux ^ext = ^DC + ^AC smiujt) threading the SQUID 
loop {^DC,AC = fji-oS H Dc^AC^ is the external flux amplitude, with //q being the permeability of 
the vacuum and S the loop area of the SQUID). The supercurrent Inm produces a magnetic field 
which couples that SQUID with its first neighbors in the x and y directions, due to magnetic 
interactions through their mutual inductances Mx and My, respectively. The dynamic equations 
for the (normalized) fluxes fnm can be written in the form |22j 

+ T^^T" + /"»^ + Sin(27r/nm) — ^x{fn-l,m + fn+l,m) — Ay(/n,m-l + fn,m+l) 

= [1 - 2{Xx + Xy)]fext, (1) 

where the following relations have been used 

T = UJot, UJo = l/VTC, fnm = ^nm/^0, fext = ^ext/^O, P = Pl/2tT = LIJ^q. (2) 

In the earlier equation, $0 is the flux quantum, Pl is the SQUID parameter, 7 is the dissipation 
constant, and Xx,y are the coupling coefficients in the x and y directions, defined as Xx^y = Mx^y/L, 
respectively. The time derivative of fnm corresponds to the voltage Vnm across the JJ of the nm— th 
rf SQUID, i.e., Vnm = dfnm/dr. The normalized external fiux f^xt is given by 

fext = foe + fAC cos (Or), (3) 

where fAc = ^Ac/^Oi foe = ^Dc/^o, and O = uj/ujo, with $dc being a constant (DC) flux 
resulting from the time-independent component of the magnetic field H. 
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Figure 2: (a) Frequency band as a function of k for a ID rf SQUID array, for /3 — 1.27, and A = —0.05 
(narrowest band, black-solid curve), A — —0.1 (red-dashed curve), A = —0.3 (widest band, green-dotted 
curve), (b) Group velocity Vg (black-solid curve) and phase velocity Vph (red-dotted curve), for a ID rf 
SQUID array with (3 = 1.27 and A = -0.1. 

The dispersion for small amplitude flux waves is obtained by the substitution of f = A exp[i(Ka;n-|- 
KyiTi — r^r)], into the linearized Eqs. ([T]) for 7 = and fext = 0, which gives 



= y I + Pl -'2{Xx cos + Xy cos Ky), (4) 

where k = {Kx,Ky) = {dx kx,dy ky). The corresponding one-dimensional (ID) SQUID array is 
obtained by setting Xy = 0, Xx = X, and by dropping the subscript m in Eqs. ([1]). 

Typical dispersion curves 17 (k) for the ID system are shown in Fig. 2a for three different values of 
the coupling A. The bandwidth AO, = ^Imax — ^min decreases with decreasing A which leads, for 
A << 1 [24J, to a nearly flat band with Ail ~ 2AV1 + Pl (and relative bandwidth ~ 2A). 

Importantly, the group velocity Vg, which defines the direction of power flow, is in a direction 
opposite to the phase velocity Vph, as it is observed in Fig. 2b. 




Figure 3: Time evolution of dissipative breathers during one period, for A = —0.1, Tb — 6.6, 7 — 0.001, 
P — 1.27, and (right panel) foe — 0.5, fAC = 0.2 - low-amplitude breather; (left panel) foe = 0, Jac = 0.6 
- high-amplitude breather. Only part of the array (iV — 30) is shown for clarity. 

3. DISSIPATIVE BREATHERS AND MAGNETIC RESPONSE 

For the generation of DDBs in rf SQUID arrays, we use the algorithm developed by Marin et. al. 
[7]. With that algorithm, we can construct low- and high-amplitude DDBs up to some maximum 
value of the coupling, Xmax, which generally depends on the external flux amplitudes /ac and 
foe [22]. Both the central site and the background of those DDBs are oscillating with frequency 
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Figure 4: Left panel: A snapshot of a two-dimensional dissipative discrete breather (DDB) for Xx = Xy = 
—0.1 and the other parameters as in the left panel of Fig. 3. Right panel: Temporal evolution of Pin 
(red-solid curve), f^'^ (green-dashed curve), and fext (black-dotted curve) during one period Tb, for (a) the 
central site of the DDB shown in the left panel {n = m = rib = N/2); (b) the site with n = m = 7 of the 
DDB shown in the left panel. 



= 2'k/Ti, = ri, i.e., the same as that of the external flux, $7. Typical single-site bright DDBs of 
both low- and high-amplitude are shown in Fig. 3 (right and left panels, respectively), where the 
spatio-temporal evolution of the induced currents i„ (n = 1, 2, 3, A^) are shown during one DDB 
period Tf,. We should note the non-sinusoidal time-dependence of the oscillations in both panels of 
Fig. 3. The linear stability of DDBs is addressed through the eigenvalues of the Floquet matrix 
(Floquet multipliers). A DDB is linearly stable when all its Floquet multipliers rrii, i = 1,...,2N, 
lie on a circle of radius Re = exp(— 7T{,/2) in the complex plane. The DDBs shown in Fig. 3 are 
indeed linearly stable. Moreover, those DDBs were let to evolve for large time intervals (i.e., more 
than 10^ Th) without any observable change in their shapes. With the same algorithm, we can also 
construct 2D dissipative breathers. A snapshot of such a DDB taken at maximum amplitude of 
the central site is shown in the left panel of Fig. 4. 

The normalized flux through the nm— th SQUID can be casted in the form 

l3inm = fnm ~ fext t (^) 

where 

fnm ~ fnm — Ax(/n-l,m + fn+l,m) — ^y{fn,m-l + fn,m+l), fext — ~ 2{Xx + Xy)]fext- (6) 

After division by the area of the unit cell cf' of the 2D array, the terms /g//, /^m> ^^'^ P^nm in Eq- 
([5]) can be interpreted as the effective external field, the local magnetic induction at the nm— th 
cell, and the magnetic response at the TiTfi — th cell, respectively. The temporal evolution of Pinm^ 
fnm: ^^'^ the external field fext, are shown in the right panel of Fig. 4, for two different sites of 
the 2D DDB shown in the left panel of Fig. 4: the central DDB site at n = m = rib = -^/2, and 
the site located at n = m = 7 (Figs, (a) and (b) of the right panel of Fig. 4, respectively). We 
observe that in the cell corresponding to the central DDB site the magnetic response is in phase 
with the applied field providing a strong paramagnetic response, while in the cell corresponding 
to the site located in the background the magnetic response is in anti-phase with the applied field 
providing moderate diamagnetic response. Thus, the local magnetic induction is sharply peaked at 
the central DDB site, as can be inferred by comparing the green-dashed curves in (a) and (b) in 
the right panel of Fig. 4. 

4. CONCLUSION 

In conclusion, we have shown using standard numerical methods that periodic rf SQUID arrays in 
an alternating external flux support low- and high-amplitude linearly stable DDBs. Those DDBs 
are not destroyed by increasing the dimensionality from one to two. Thus, we have constructed 
several linearly stable DDB excitations both for ID and 2D rf SQUID arrays, which may alter 
locally the magnetic response of the arrays. Planar SQUID arrays similar to those described here 
have been actually constructed and studied with respect to the ground state ordering of their 
magnetic moments [23]. Thus, the above theoretical predictions are experimentally testable. 
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